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An experimentally validated approach is described for fast axisymmetric 
Stirling engine simulations. These simulations include the entire displacer 
interior and demonstrate it is possible to model a complete engine cycle in 
less than an hour. The focus of this effort was to demonstrate it is possible 
to produce useful Stirling engine performance results in a time-frame short 
enough to impact design decisions. The combination of utilizing the lat- 
est 64-bit Opteron computer processors, fiber-optical Myrinet communica- 
tions, dynamic meshing, and across zone partitioning has enabled solution 
times at least 240 times faster than previous attempts at simulating the 
axisymmetric Stirling engine. A comparison of the multidimensional re- 
sults, calibrated one-dimensional results, and known experimental results 
is shown. 

This preliminary comparison demonstrates that axisymmetric simula- 
tions can be very accurate, but more work remains to improve the simula- 
tions through such means as modifying the thermal equilibrium regenerator 
models, adding fluid-structure interactions, including radiation effects, and 
incorporating mechanodynamics. 
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1. Oscillating flow which changes the effective conduction, flow friction coefficients and 
heat diffusion. 68 

2. Low mach number flow (no shocks), 

3. Compressible flow due to enclosed varying volumes and heat transfer effects, 

4. Laminar, Transitional, and Turbulent flow with Reynolds numbers from 100 to 10000 
(based on various length scales pertinent to the flow region), 

5. Conjugate heat transfer, thermal dispersion and local thermal nonequilibrium in the 
porous media regenerator, 

6. Micron to Millimeter scale geometry 6 

7. Sliding Interfaces in the appendix gap 

8. Deforming Flow Regions Due to Compression and Expansion 


III. Current Multi-Dimensional Stirling Cycle Analysis Tools 


At the level of multidimensional analy- 
sis, relatively little has been completed be- 
cause the third order analysis, including 
Sage, 5 and to some extent, DeltaE, 13 and 
MS*2 Stirling Cycle Code 11 are faster and 
for the most part have been adequate engi- 
neering design tools. However, to improve 
efficiency, reliability, and supplement exper- 
imental testing (to understand and reduce 
losses, identify component temperatures, as- 
sist in sensor calibration, correlation, and 
placement) it will likely require a better un- 
derstanding of the actual flow features and 

heat transfer throughout the engine. Some of the primary multidimensional analysis tools 
are shown below. 



Figure 3. Example Transient Simulation of 
Stirling Convertor-Colored by Temperature 


1. Modified, Computer Aided Simulation of Turbulence (CAST) 

The modified CAST code 70 is based upon the Semi-Implicit Method for Pressure-Linked 
Equations (SIMPLE) 71 method but is restricted to two dimensions. 72-74 It was modified to 
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Nomenclature 


7 Porosity of the Medium 

r e ff Stoke’s Tensor 

pf Fluid Density 

p s Solid Density 

v Fluid Velocity Vector 

Ef Fluid Energy 

E s Solid Energy 

h Sensible Enthalpy 

J Diffusion Flux of Helium 

k e ff Effective Thermoconductivity 

p Pressure 

Sf Fluid Enthalpy Source Term 

T Temperature 


I. Introduction 

Great strides have been made in the last decade in expanding the designer’s analy- 
sis options for building improved free-piston Stirling engines 1 ranging from thermal and 
structural dynamic analysis, 2,3 computer-assisted-design (CAD) software for complete three- 
dimensional part design, 4 one-dimensional thermodynamic cycle analysis 5-18 and optimiza- 
tion, and one-dimensional system mechanodynamic simulations. 19-25 A review of the progress 
in Stirling analysis shows a steady relaxation of assumptions/limitations imposed on the 
analysis over the past century. 26 

Initial two-dimensional computational fluid dynamic (CFD) simulation work began over 
twenty years ago and focused on specific regions of the engine but were not time-efficient 
and the projects were not continued. 6,27-30 And other researchers attempted to build specific 
in-house research codes with mixed success at modeling specific components. 31-36 

More recently, multidimensional multiphysics simulations that allow for realistic Stirling 
engine system solutions now exist and include the following capabilities that are potentially 
useful for Stirling convertor design: 37-40 

• Structural 

— linear or nonlinear static or transient deformation 

— modal, spectrum, harmonic, and random vibration dynamics 

— linear or nonlinear buckling 
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• Thermal 


— steady-state, transient, conduction, convection, and radiation 


• CFD 


— Steady-state/transient, compressible/incompressible, laminar /turbulent, forced/natural 
convection, conjugate and radiation heat transfer, laminar to turbulent transition 
modeling, scalable parallel performance 

• Electromagnetics 

— electrostatics, magnetostatics, low-frequency electromagnetics, current conduc- 
tion, circuit analysis, harmonic, transient, modal, scattering, perfect electric and 
magnetic conductors, impedence boundaries, near and far electromagnetic held 
extension, frequency selective surface, antenna radiation patterns, and specific 
absorption rate. 

• Coupled Physics 

— acoustics-structural, electric-magnetic, fluid-structural, fluid-thermal, electromagnetic- 
thermal, piezoelectric, piezoresitive, thermal-electric, thermal-structural, electromagnetic- 
thermal-structural, and eletrostatic-structural 

Recently, attempts at utilizing these commercially available options have focused on only 
the thermal dynamic cycle simulation because it is difficult to achieve a useful solution 
in a reasonable time, let alone attempt a more comprehensive simulation. For example, 
previous attempts at whole engine simulation with these multiphysics commercial codes 
ended in unconverged solutions. 31,41 Although a recent simplified axisymmetric solution 76 
shows reasonable energy balance convergence, the results are not validated yet by experiment. 

The first reported successful three-dimensional whole Stirling engine simulation was re- 
ported by Mahkamov 42-45 in the United Kingdom. Two gamma-type Stirling engines were 
simulated using commercially available software with good reported experimental agreement 
of 12%-18% error in the indicated power. Seal leakage, appendix gap flow, and other loss 
mechanisms were not included in this work. 

The second reported successful three-dimensional Stirling engine simulation was by Zhang 46 
in China in which it took approximately 3 months to reach a steady-harmonic converged 
solution on a simplified, academic geometry that did not include the appendix gap, solid 
walls, or internal displacer. 

A multi-dimensional analysis has an important place in Stirling engine design as is dis- 
cussed at length in Ref. 47 and outlined below: 
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• verify the one-dimensional results, 

• properly simulate inherently three-dimensional turbulence-including transition, 

• provide empirical heat transfer and friction factor coefficients for complex geometries 
before the hardware is built 

• integrate all the parts into a CAD system and test for structural and relative motion 
issues in the overall design 

• assist experiments with determining hard to reach data 

• provide a fluid-structure interaction capability 

• generate linear reduced order models for controller development 

• solve high-power Stirling applications in which the one- dimensional flow assumption 
breaks down 

• identify areas of excessive flow losses due to unintended dead zones, recirculation zones, 
dissipative turbulence and other losses such as: 48-50 

1. Inefficient heat exchange and pressure loss in the regenerator, heater and cooler, 

2. Gas spring and working space loss due to hysterisis and turbulence, 

3. Appendix gap losses due to pumping and shuttle effects, 

4. Mixing gas losses from nonuniform temperature and flow distributions perpendic- 
ular to primary engine flow axis, 

5. Conduction losses from the hot to cold regions 

6. Losses due to combined radiation, conduction and convection in hollow spaces 

7. And in general, inaccurate loss representations due to use of 1-D flow design codes 
to account for flow and heat transfer through area changes (between components) 
where phenomena such as flow separation and jetting from tube into regenerator 
may occur. 

The current multiphysics simulation tools are not optimized for Stirling analysis as dis- 
cussed in Ref. 51 For example, recently developed oscillating time advance strategies 52-58 and 
high order spatial differencing could make modeling oscillating turbulent transition a reality. 
But given the cost of developing new codes, it is desirable to fully utilize current technology 
first. 
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This paper will present the first successful U.S. axisymmetric Stirling engine simulation 
incorporating an actual engine’s geometry including appendix gap, displacer seals (inner and 
outer), solid walls, and displacer interior. Only the flexures, piston seals, bounce space, radi- 
ation shields, and heat exchanger fins were not included for this initial work, but reasonable 
agreement with experiment has been achieved despite these omissions. 

II. Description of the Problem 

The dual opposed configuration shown in figure 1 59-62 is being developed for multimission 
use (i.e., for use in atmospheres and space), including providing electric power for potential 
missions such as unmanned Mars rovers and deep space missions. 63,64 

Only the Stirling engine part of the 
convertor (Fig. 2) is simulated multi- 
dimensionally although one could anticipate 
the entire convertor may one day be simu- 
lated using modified multiphysics software. 

An example of the desired full Stirling 
engine simulation is shown in Fig. 3. It is 
anticipated that full 3D simulations will pro- 
vide a level of geometric and flow detail necessary for further design improvements. For 
example, hardware experiments have shown that large performance gains can be made by 
varying manifolds and heat exchanger designs to improve flow distributions in the heat ex- 
changers. 60,65-67 



Figure 1. Dual Opposed Stirling Convertors 
Reduce Vibration (Schreiber) 



Figure 2. Stirling Engine Part of Convertor Simulated Only 

The kinds of flow and geometry that occur in a Stirling engine for which a modeling 
technique must address are as follows: 49,50 
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include oscillatory boundary conditions and conjugate heat transfer but does not allow for 
complicated geometry nor sliding interfaces. It has been used to model Stirling components 
but has not been extended to a whole engine simulation tool. 

A pressure-splitting technique was added 75 to reduce the computational requirements. 
It was based on separating the thermodynamic and hydrodynamic pressures so that these 
widely varying scales could be solved with less round-off error and better efficiency. And the 
non-acoustic form of the equations was later incorporated into the code since the fast acoustic 
waves and small dimensions of the engine essentially resulted in ’’incompressible” like flow 
behaviour and removing the acoustical pressure disturbance could enhance computational 
speed and accuracy . 32 

2. CFD-ACE 

This commercial code has been used to model a two-dimensional representative Stirling 
engine 31, 33,46,76 and it has been used to model an actual engine . 41 Full validation of these 
simulations has not been completed. 

It is also based upon the SIMPLE technique. The regenerator is not currently modeled 
correctly since thermal equilibrium is assumed between the gas and solid, but as shown 
later, this error may not be large compared to the overall engine performance. The code 
does support sliding interfaces on parallel computers in three-dimensions and the company 
is introducing sliding, overlapping Chimera grids which should enable better appendix gap 
modeling. Moreover, the highly compressed volumes occuring at both ends of the displacer 
may be more easily modeled with Chimera grids. This finite volume code can utilize both 
structured and unstructured grids. 

3. Fluent 

This commercial code is based upon the SIMPLE and PISO methods. It currently has similar 
regenerator modeling limitations as CFD-ACE. It does however have a sliding interface for 
two-dimensional axisymmetric flows that can be used for appendix gap modeling on parallel 
computers and will be discussed in more detail later. It is also finite volume based and can 
utilize both structured and unstructured grids. 

It has been used with mixed success by industry for Stirling engine simulations 77 by 
several commercial manufacturers. It was also used by Mahkamov 43 in what appears to be 
the first full three-dimensional simulation, but details are embargoed. Also, in the related 
business of cryocooler modeling it has been used with reasonable success on simplified prob- 
lems 78-80, 82-89 And this report has utilized it to provide a detailed axisymmetric simulation. 
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I STAR- CD 


The Simulation of Turbulent Flow in Arbitrary Regions (STAR) 38 code also uses SIMPLE 
and PISO methods. Its companion product (STAR-HPC) is the parallel computer version. It 
also has sliding interfaces and deforming mesh capability. It has been used in the related held 
of internal combustion engine piston modeling, and some Stirling engines and regenerators 
have been modeled with it. 90,91 

5. CFX 

This code also uses SIMPLE and PISO methods on unstructured grids. It also has sliding 
interfaces implemented, but no Stirling engine modeling with this software has been publicly 
published. 37 It currently seems to have the greatest potential for full multiphysics Stirling 
convertor simulations because of it’s fluid-structure interaction, magnetic flux analysis, and 
structural/thermal stress analysis capabilities. It’s dynamic meshing capabilities are not 
ideal for highly compressed volumes as occur in the Stirling engine since only the spring 
analogy mesh adaption is currently available. 

6. Others 

While there are other research codes, they are usually limited to modeling only specific 
regions of the Stirling engine such as the regenerator, or the displacer. 15, 81,82,92,93 

IV. Regenerator Modeling 

A very important specific area of modeling difficulty is the regenerator. Since the regen- 
erator (depending upon one’s definition of effective) has roughly 3 to 40 times more effective 
heat transfer than the heater, 94 any inefficiency of the regenerator represents a significant 
loss for the entire Stirling engine. Hence, any numerical losses or inaccuracies in this region 
will disproportionately influence the entire Stirling simulation. 

Also, the geometrical shape of the matrix elements in the regenerator is important and as 
shown in Figs. 4(b-f), many shapes are possible. For practical calculations a Representative 
Elementary Volume (REV) is often used as shown in Fig. 4(a). An REV requires some 
closure assumptions that can be partially derived from experiment and potentially from 
microscale (pore level) numerical simulations. 

Since the size and shape of the regenerator fibers do affect Stirling engine performance, 
numerically simulating a regenerator has been pursued at both microscopic (direct) and 
macroscopic (REV) levels of fidelity. 6, 95-101 An example simulation utilizing Fluent by Har- 
vey from Georgia Tech is shown in Fig. 5. 
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(a) Representative Elementary (b) Random (c) Wire Mesh 

Volume 



(d) Perforated Disk (e) Foam Metal (f) Sintered Glass (Failed) 


Figure 4. Regenerator Geometry 


However, our initial approach at whole engine simulation simply adds source terms to the 
governing equations and does not take into account the different temperatures of the solid 
matrix and the fluid. The source terms serve to enforce a pressure drop according to the 
Darcy-Forcheimer 68 equation. And only a single energy equation is used, although it does 
include an effective solid/fluid average conductivity 69 and an effective solid/fluid average 
energy. 

The continuity equation is not modified, but the following source term is added to the 
momentum equations (for i — x,y, orz ): 


Si = 



V T f- 2 2 P^rnag ^i 


(1) 


where permeability K = 4.08e -lo m 2 , C2 = and inertial coefficient, Cf = 0.178. 
These coefficients were determined experimentally 102 and used without modification here. 
And the fluid/solid averaged energy equation can be written as: 
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(a) Schematic (b) Pressure (c) Pathlines 

Figure 5. Flow Through Wire Mesh Screen Regenerator 


V. Dangers of Component Modeling 

Modeling the entire engine, as demonstrated herein, is preferable to modeling one com- 
ponent at a time due to the oscillatory nature of the flow. A component requires known 
unsteady boundary conditions for an entire cycle. Since a given boundary may have flow 
going in both directions, the direction of the characteristics at these boundaries are diffi- 
cult to define. Incorrect specifications of boundary conditions result in solution instability, 
nonconvergence of solutions, and/or convergence to inaccurate results. It is well known 
that even for steady harmonic flows, non-physical reflections can occur at inflow and out- 
flow locations. 103 In essence, by artificially truncating the domain to model a single engine 
component an artificial boundary condition is required. However, the exact specification of 
artificial/non-reflecting boundary conditions in nonlinear turbulent flow is an important and 
as yet unresolved area of research. 104 

For example, a subsonic inflow has four characteristics entering and hence four physical 
boundary conditions need to be specified, typically total pressure, total temperature, and two 
flow angles, allowing only one of the velocity components to be free (which is extrapolated 
from the interior flow solution) . A subsonic outflow has only one characteristic entering and 
hence only one physical boundary condition should be specified, for example, static pressure. 

Simultaneous inflow/outflow along each engine component’s interface is clearly demon- 
strated by examining flow velocities near zero (from -1 m/s to 1 ms). Regions of blue and 
red in Fig. 6 show borders of mixed inflow/outflow in the axial direction. For example, the 
piston compression space in Fig. 6(a), the seal exit in Fig. 6(b), the cooler exit in Fig. 6(c), 
and the appendix gap entrance in Fig. refflg:oscillating(d) all show mixed inflow/outflow 
borders at this moment in time. Moreover, these mixed boundaries moving during the cycle, 
further complicating the numerical boundary treatment. If a component boundary crosses 
any of these indicated areas, the boundary conditions will not be properly defined. Of course 
in the case of commercial software users the code will attempt to stabilize the solution with 
essentially guessed boundary conditions, but the solutions will be inherently non-physical. 
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(a) Entire Engine 


(b) Seal Region 




(c) Compression Space 


(d) Heater Manifold and Appendix Gap 


Figure 6. Oscillating Axial- Direct ion 


Repeating the same velocity test, but this time examining velocity close to zero in the 
y-direction, we see a new set of boundary lines that will pose difficulties resulting in non- 
physical component solutions in Fig. 7. 


A. Adjust Components within a Full Engine Simulation 

Instead of pursuing the mathematically ill-posed problem of component modeling in com- 
plicated oscillating flow, it is easier, more efficient, and more reliable to determine a whole 
engine steady-harmonic solution with the procedures described herein. The engine compo- 
nents can be easily examined and it only takes a few hours to test the effect of component 
changes since the overall engine has the correct steady-harmonic temperature distribution. 

At the same time, it is possible to numerically determine the one- dimensional equivalent 
friction factors and heat transfer terms in the new components required by one-dimensional 
analysis tools (e.g. Sage). In this way, a one-dimensional analysis be made more accu- 
rate while not requiring the ill-posed boundary conditions that multidimensional component 
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(a) Entire Engine 


(b) Heater Region 



(c) Appendix Gap Region (d) Seal Region 

Figure 7. Oscillating Radial-Direction 

modeling entails . 47 Finally, it is worth noting that one-dimensional analysis does not have 
difficulty with component modeling because a boundary is always either an inflow or an 
outflow, never both. 
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VI. Approach for Fast Simulations 


Achieving practical computer simulations of entire Stirling engines depends upon judicial 
use of computer hardware, software engineering, grid generation, and algorithm design. A 
combination of all four areas is presented next that enables axisymmetric whole engine 
simulations to be accurately completed in less than a week. 


A. Computer Design 


An important area new to Stirling analysis is the use of parallel computers for faster simu- 
lations. NASA Glenn was one of the first organizations to pioneer computer clusters back in 
1992. Many universities have assembled their own low cost computers, 105 and many states 
have supercomputer consortiums 106 to assist in designing and operating small computer clus- 
ters that utilize commercial off-the-shelf (COTS) hardware and software for low cost systems. 
These systems are referred to as Linux Beowulf clusters 107, 108 and now represent the ma- 
jority of the top 500 supercomputer systems in the world. 109 Recently, many commercial 
organizations have begun assembling ’’Beowulf’ -like clusters including Dell, 110 Angstrom, 111 
LinuxNetworx, 112 Appro, 113 and Microway 114 as shown in Fig. 8. They cost a bit more, but 
provide comprehensive health monitoring and controls that significantly reduce the admin- 
stration and integration costs. The latest systems can house 260 processors in a single 7 foot 
tall tower. 



(a) Angstrom Cluster 



Space Tar 
2 switches 


Power reed 


Rear View 


1 master node 


(b) Appro Cluster 



(c) Linux Networx Cluster 


Figure 8. Appro, Angstrom, and Linux Networx Clusters 


Our sixteen dual-Opteron (32) processor Microway cluster is shown in Fig. 9. This was 
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selected because last year it represented the best value for performance and it integrated 
well with the commercial software, Fluent. The Opteron chip, at the time, was the only chip 
that could run 32-bit and 64-bit applications with an exceptional price/performance ratio. 
All the whole Stirling engine results shown here were calculated using it. 



(a) Front View (b) Dual CPU Board 

Figure 9. Microway Cluster 


A very important consideration when designing your computer cluster is the commu- 
nication bandwidth and latency since commercial fluid simulation software utilizes older 
numerical approaches that require a large number of frequent communications between the 
compute nodes. In Fig. 10, two of the best communications network interface card (NIC) and 
network switches are shown. The Infiniband system can perform 10 Gb/s bi-directional com- 
munication rates or about 20 times faster than the Myrinet in our cluster. Our simulation 
utilized Myrinet communications because last year Infiniband was not a proven technology. 



(a) Myrinet Card 


(b) Myrinet Switch 


(c) Infiniband Card 


(d) Infiniband 144 port 
Switch 


Figure 10. Communications 
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B. Software Decision 


We chose Fluent as our computational platform because, at the time, it was the only com- 
mercial vendor supplying an axisymmetric sliding interface for parallel computers. Presently, 
as mentioned in Sec. Ill, the other vendors are now also providing similar capability. It is 
important to always have at least two commercial solvers available because of the invariably 
large number of bugs in their software. Generally, a bug is more easily found if the same 
simulation is repeated with two different codes. For example, Fluent version 6.1.22 properly 
handles heated moving solids, but the new release (version 6.2.16) did not as identified by 
Ibrahim. 115 The latest version 6.2.19 has been corrected. 

Another reason for utilizing at least two commercial codes is the lack of control the user 
has in the solvers utilized. Generally, the commercial codes all utilize slight variants of 
older, well established lower order numerical techniques. The slight variants are the only 
independent way to examine the effect of a given solver and to look for possible defects. 
Of course developing one’s own software offers more freedom, but this can be risky in the 
short-term. 

Our second platform utilizes CFD-ACE and reasonable solutions have been achieved with 
it. It’s primary limitation has been the need to work in three-dimensions when solving cases 
including axisymmetry sliding interfaces. Another strength of this software is it includes a 
wide range of multiphysics capabilities. 

A third very interesting platform is CFX because of it’s integrated multiphysics capability 
enhanced by the merger between ANSYS and CFX. The only apparent limitation of the 
software is grid layering is not yet available. Grid layering is defined by the automatic 
addition or removal of a row of cells as the mesh domain shrinks or expands and will be 
discussed later. 

Other packages such as STAR-HPC are very capable as well. However, the commercial 
cost of these packages are very expensive and it is not practical to test all of them. 

Finally, full three-dimensional simulations are still quite challenging and the parallal 
scalability of the commercial software is somewhat limited due to the implicit nature of 
their steady flow solvers. Better approaches are currently available 51 but have not yet been 
applied to the Stirling cycle analysis. 

C. Grid Design 

When using commercial codes, it is primarily the grid quality that determines the speed 
and accuracy of the solution. Ideally, one could simply read in a three-dimensional CAD 
geometry and a grid would automatically be created. That is in fact the implied goal of 
computational fluid dynamics code developers. But for complicated geometry as occurs 
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inside a Stirling engine with all it’s moving parts, the technology currently available still 
requires a great deal of human intervention. Taking several weeks to develop a grid is not 
unusual and is often required to get a reasonable solution. For example, the rate at which 
conjugate heat transfer can be solved is directly related to the smoothness and quality of 
the grid at the interface. This is particularly true in Stirling engines when the heat capacity 
of the solid is vastly diffent than the Helium. If short-cuts are taken in developing the grid, 
the price is often vastly longer solution times, if a solution is produced at all. 

All the grids produced in this work were done with Gambit, which is the Fluent, Inc. 
proprietary grid generator. There are other grid generators available with their own distinct 
advantages, but despite the bugs encountered in Gambit, we found it to be a useful tool. 

First, the domain is scaled up by 1000 to compensate for the 32-bit precision available in 
the grid generator. The Stirling engine, when entirely modeled, has a fairly wide range 
of length scales ranging from a thousandth of an inch in the narrow seal and appendix 
gap region to an inch or more in the heat exchanger regions. This scaling has the effect 
of making very small numbers larger so that smoother grids can be created. Once the 
mesh is created satisfactorly, the entire grid domain is scaled back down by a 1000 
before utilizing it. 

Second, try to use structured or paver unstructured grids wherever possible to reduce the 
number of bad cells and to reduce the number of cells required. Avoid unstructured 
tetrahedrals because of the considerably reduced rate of solution convergence. 

Third, be sure a sliding interface has essentially matching cell sizes on both sides of the 
interface to prevent interpolation errors and avoid instabilities. 

Fourth, be sure to place grid points in the boundary layer to properly model friction losses. 

Fifth, use sizing functions to smoothly expand the cell size from the boundaries at a rate 
of 10 percent. 

Sixth, use layering wherever possible, but avoid grid adaption and remeshing because this 
slows down the solution time considerably. Choose a layering size large enough to 
permit a reasonable 160 time steps per cycle, but small enough to allow for smooth 
transitioning to boundary layers. 

Seventh, utilize moving boundary layers on all walls that are moving to enhance conjugate 
heat transfer. And place a boundary layer on both sides of the interface (on the solid 
side and the fluid side). 
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Eighth, include as much geometry as is practical to avoid having to redo the grid later when 
simplifying assumptions are found to matter. 

Ninth, clean up the geometry first in a CAD package to avoid creating virtual geometry 
that is not as robust. 

Tenth, test your grid by forcing the piston and displacer to be stationary and simulate the 
heat transfer in the engine. Once a converged solution is observed, continue running 
the solver and look for grid induced heat sources. It should be noted that Cleveland 
State University found producing a steady-state result first provided a better starting 
condition for the transient solution. But this last suggestion goes a bit further in 
providing a means for improving the grid by ’’over-converging” the steady solution. 


For the successful case presented here, we used about 700,000 grid points and all grid 
boundary layers start at a distance of .001 in from the fluid/solid interface. 




(a) Entire Engine 


(b) Piston Region 



(c) Spider Space 


Figure 11. Grid Strategy 
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In Fig. 11, we show the entire engine, piston region and spider region grid utilized in 
this simulation. Notice the piston has a moving boundary layer and grid layering is used 
to avoid small cells when the piston fully compresses. Also, notice in the spider region 
that all helium/wall interfaces have boundary layers on both sides to enhance conjugate 
heat transfer. The overall engine can be seen to not include any grid inside the heater or 
cooler since the heat exchanger wall temperture is specified as constant. Also notice there 
is some extra wasted grid in the spider region near the compression space. This is necessary 
to maintain grid quality in the other nearby boundary layer regions. Grid quality is more 
imporant than grid count for faster conjugate heat transfer calculations. 

In Fig. 12, the hot half of the engine is shown, including the gridding in the displacer and 
appendix gap. Notice that despite the complexity, tetrahedral grids were not required. The 
layering inside the displacer was done to avoid moving all the grid points inside the displacer 
for faster run-times. For example, without grid layering (the addition or subtraction of rows 
of cells), all the grid cells inside the displacer would need to move. Instead, only the single 
layer of cells located at the top of the displacer (right side in picture) moves at a time in 
Fig. 12(b). Also, the expansion space uses layering to avoid skewed cells when the displacer 
is at top dead center (TDC). 

In Fig. 13, the gridding utilized in the seals and compression space are shown. The 
challenge in the seal region is to avoid overly skewed cells, overly mismatched cells across the 
interface, and reasonable cell size variation near the boundary layer. It is quite difficult to 
balance all these constraints. Layering was used in the compression region to avoid crushed 
cells and the axial direction spacing was chosen to maximize the time-step. 

In Fig. 14 the seal inside the displacer is shown. Notice the use of unstructured quadri- 
laterals throughout to avoid the inefficient tetrahedral gridding. Also notice just how small 
the seal gap is even compared to the thickness of the displacer wall. This region was very 
difficult to grid and required utilizing suggestion 10 in Section C. Velocities are high in this 
region and this grid may not be adequate here, but since this seal region is not being studied 
this was left as seen. 
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(a) Right Half 


(b) Layering Inside Dome 




(c) Layering in Expansion Space 


(d) Appendix Gap 



(e) Appendix Gap Close 


Figure 12. Hot End Grid Strategy 
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(a) Top and Bottom Seal, Rod, Compression Space (b) Top Seal and Compression 



(c) Top Seal Close-up (d) Bottom Seal and Compression 

Figure 13. Seal Grid Strategy 
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(c) Inside Seal Close-up (d) Inside Seal Close-up2 

Figure 14. Inner Seal Grid Strategy 
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D. Solver Optimization 

While one’s choices are limited with commercial codes, they do offer a variety of low-order 
solvers to choose from, including low Mach number variable density flows. 116 A comparison 
of two commercial codes with the latest high-order numerical techniques demonstrates some 
of the differences a solver decision can make. 51 One of the limitations currently in Fluent 
is only first order accuracy is available in time when the mesh moves. CFD-ACE permits 
a higher accuracy in time with moving meshes, but does not have a coupled solver. The 
results shown in the next section utilize the segregrated implicit solver in Fluent, although it 
was later discovered the coupled implicit solver is potentially significantly better. Of course, 
the grid used for the segregrated solver may need further improvement if a coupled solver 
is used. It was found that the best method in Fluent and the best method in CFD-ACE 
performed comparably on simpler heat transfer test problems. 51 And with sufficient grid 
quality and density, the answers are comparable to the latest numerical methods. 

One of the challenges to utilizing actual engine geometry is the difficulty of creating a 
high quality mesh as discussed previously. However, certain switches in the solvers may 
be adjusted to attempt to compensate for poor grid quality in regions of conjugate heat 
transfer. 117 In what follows are specific recommendations for improving conjugate heat 
transfer and they contain Fluent specific switches. The other commercial packages should 
be modified with their own specific format. 

First, skewed meshes require secondary temperature gradients to be turned off to improve 
stability and to reduce undershoot/overshoot of local cell temperatures. In Fluent, this 
is accomplished by the command, ’’(rpsetvar ‘temperature/secondary-gradient? #f)” 

Second, in conjugate heat transfer problems with vastly different thermal conductivities 
between the fluid and solid zones, cells near the interface can experience round-off 
errors. To mitigate this effect, we want to ensure that the coarse levels of the AMG 
are visited. It is recommended that the default AMG cycle for energy be changed from 
Flexible to either F, V, or W-cycle. The F and V-cycles are more efficient in parallel 
simulations. The termination criteria should also be reduced to either 0.01 or 0.001 in 
order to tighten up convergence. 

Third, although a reasonable temperature distribution can be achieved using single pre- 
cision, in order to ensure that the total heat transfer is balanced the use of double 
precision is recommended. 

Fourth, if convergence still remains problematic with secondary gradients turned off, it is 
recommended that explicit under-relaxation of temperature be used to stabilize the 
solution. When this is activated, the energy equation is solved with a URF of 1, but 
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the solution is explicitly relaxed with respect to the previous iteration. This allows 
for the equations to be efficiently solvered (URF equal 1) and instabilities due to poor 
quality mesh and non-linear material properties to be mitigated (relaxation factor). 
The GUI entry of URF for energy will specify the explicit relaxation factor to be used. 
To active this, the scheme command is: (rpsetvar ‘temperature/explicit- relax? #t) 
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VII. Comparison with Analytical and Experimental Data 


The approach just described was applied to simulating a Stirling engine built by Infinia 
Corp. 118 (formerly Stirling Technology Company) being considered for use in space missions. 
The engine is currently being tested for long-duration missions and hence a good deal of 
experimental performance information is available for comparison. The geometry is export 
controlled for these engines and so limited information may be shared with the public. The 
entire simulation was performed first, with no knowledge of the ’’correct” answer either from 
experiment or calibrated one-dimensional analysis. The operating conditions used in the 
simulation are shown below (Table 1):. 


Hot-End Temperature, K 

923 

Cold-End Temperature, K 

353 

Ambient Temperature, K 

293 

Frequency, Hz 

80 

Mean Pressure, Pa 

2.429E+6 

Power Piston Amplitude, mm 

6.0e-3 


Table 1. Operating Conditions 


The pressures and heat transfer occuring in various parts of the engine over two cycles 
are shown in Fig. 15(a). The pressure in the left (cooler end), middle, and right (heater 
end) of the regenerator is shown in Fig. 15(b). The average pressure in the expansion and 
compression spaces are shown to be about equal in Fig. 15(c). And, the temperature after 
412.376 cycles of simulation are shown in Fig. 16 for an instant in time (135.36 degree Crank 
angle with reference to 0 degree displacer located at mean position) in the cold end of the 
machine and for the hot end are shown in Fig. 17. Notice the effects of radial heat transfer 
on the regenerator. 

An initial comparison with Sage suggested the axisymmetric simulation was over-predicting 
the indicated power as shown in Table 2. And the overall heat transfer in and out fell within 
a band predicted by Sage in which heat loss to the environment was included and then 
assumed excluded (adiabatic). 

But then we gathered the actual performance of two identical engines identified as Tech- 
nology Demonstrator Convertor (TDC) #13 and #14 from over many thousands of hours of 
operation since August 2003. As shown in Table 3, the axisymmetric solution which assumed 
no ambient heat loss to the environment, predicted not only the higher indicated power more 
accurately, but the efficiency was essentially within the variability of each built machine. 

A few words of caution about this encouraging result are in order. First, due to experi- 
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Absolute Pressure (Pa] 



Time 



(a) Pressure and Heat Transfer over Two Cycles (b) Pressure Along Regenerator over Two Cycles 



(c) Comparison of Expansion and Compression Pres- 
sure 

Figure 15. Pressure and Heat Transfer over Two Cycles 
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Axisymmetric 

Simulation 

Sage ID 

Ambient Heat Loss 

Sage ID 
No Heat Loss 

PV Power, W 

79.65 

68.77 

70.14 

Heat In, W 

247.315 

260.94 

193.8 

Heat Out, W 

168.313 

191.9 

123.7 

PV Efficiency 

.322 

.264 

0.362 

Pressure Ratio 

1.209 

1.187 

1.187 

Regenerator Delta Press., Pa 

10466.3 

15,888 

15790 

Heater Delta Press., Pa 

682 

194 

192 

Cooler Delta Press., Pa 

896 

316 

315 

Pressure Amplitude., Pa 

225269 

203100 

203200 

Mean Pressure, Pa 

2429130 

2378000 

2378000 


Table 2. Comparison With Sage ID Results 


mental limitations the piston amplitude is not known for certain. The original designs call 
for a piston amplitude of 6+/- .08 mm. And the mean operating pressure is 2.517 Mpa 
compared to the simulated mean pressure of 2.42913 Mpa. If the piston amplitude in the ex- 
periment is exactly 6.00mm, then this mean pressure difference implies the simulated power 
would need to be corrected up to 82.5874W. But since the engines are tuned for maximum 
performance and the energy measurements are considered more reliable than the piston am- 
plitudes, it seems appropriate to conclude this simulation did in fact predict the performance 
and efficiency within the bounds of variability of the two engines tested. Since a slightly 
lower piston amplitude would counter-balance the corrected higher power just shown, it is 
important to improve the experimentally derived piston amplitude for future comparisons. 



T h 

T 

freq 

net Qin 

PV Pout 

net Qout 

TDC #13 

646.0 

92.4 

81.4 

242.1 

78.2 

163.9 

TDC #14 

646.5 

94.4 

81.4 

250.4 

79.6 

170.8 

Simulation 

650.0 

80.0 

80.0 

247.3 

79.7 

168.3 

%Err-AVE 

0.6% 

0.3% 

-1.8% 

0.4% 

0.9% 

0.6% 

%Err-TDC#13 

0.5% 

0.2% 

-1.8% 

2.1% 

1.9% 

2.7% 

%Err-TDC#14 

0.6% 

0.3% 

-1.8% 

-1.2% 

0.0% 

-1.5% 


Table 3. Comparison With Experimental Results 
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3.75e+02 
3.74e+02 
3.73e+02 
3.71 e+02 
3.70e+02 
3.69e+02 
3.68e+02 
3.67e+02 
3.66e+02 
3.64e+02 
3.63e+02 
3.62e+02 
3.61 e+02 
3.60e+02 
3.58e+02 
3.57e+02 
3.56e+02 
3.55e+02 
3.54e+02 
3.53e+02 
3.51 e+02 
3.50e+02 
3.49e+02 
3.48e+02 
3.47e+02 
3.45e+02 
3.44e+02 
3.43e+02 
3.42e+02 
3.41 e+02 



™ 3.40e+02 

Contours of Static Temperature (k) (Time=5.1547e+00) 



(a) Left Half 


(b) Compression and Cooler 



(c) Cooler and Seal 



(d) Seal Region 



(e) Seal Pretty Close (f) Seal Very Close 


Figure 16. Cold End Temperature 
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(a) Right Half 


(b) Below Heater 




(c) Around Heater 


(d) Regenerator 



Figure 17. Hot End and Regenerator Temperature 
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VIII. Conclusion 


For the first time in the U.S., a fully converged axisymmetric simulation of an actual 
Stirling engine was successfully performed. And moreover, by utilizing the latest advances 
in computer processors, communication hardware, parallel decomposition strategy, dynamic 
meshing, solution algorithms, and careful grid design, it is possible to compute one complete 
Stirling cycle in less than an hour. 

The full axisymmetric simulation predicted the engine performance within the bounds 
of both engine test’s data. And predicted the overall efficiency and indicated power consid- 
erably better than our initial efforts with using a leading one- dimensional Stirling analysis 
prediction on the first attempt. The one-dimensional analysis can certainly be improved 
by modifying various friction and heat transfer factors, either from experimental data or 
from this axisymmetric solution. But the real value of the axisymmetric solution is no ad- 
justments were required because the engine is simulated from first principles. Of course 
the one-dimensional analysis is valuable for other reasons such as fast optimizations and 
performance mapping. 

These results and observations lead to the conclusion that one, two, and three-dimensional 
modeling should all be employed and all three paradigms provide important capabilities 
that when combined provide a potent combination of initial design, empirical coefficient 
adjustment, optimization, and final prototype demonstration before the first metal is cut. 
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